Bose-Einstein condensates under a spatially-modulated transverse confinement 
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We derive an effective nonpolynomial Schrodinger equation (NPSE) for self-repulsive or attrac- 
tive BEC in the nearly-lD cigar-shaped trap, with the transverse confining frequency periodically 
modulated along the axial direction. Besides the usual linear cigar-shaped trap, where the peri- 
odic modulation emulates the action of an optical lattice (OL), the model may be also relevant 
to toroidal traps, where an ordinary OL cannot be created. For either sign of the nonlinearity, 
extended and localized states are found, in the numerical form (using both the effective NPSE and 
the full 3D Gross-Pitaevskii equation) and by means of the variational approximation (VA). The 
latter is applied to construct ground-state solitons and predict the collapse threshold in the case 
of self-attraction. It is shown that numerical solutions provided by the one-dimensional NPSE are 
always very close to full 3D solutions, and the VA yields quite reasonable results too. The transition 
from delocalized states to gap solitons, in the first finite bandgap of the linear spectrum, is examined 
in detail, for the repulsive and attractive nonlinearities alike. 

PACS numbers: 03.75.Ss,03.75.Hh,64.75.+g 



I. INTRODUCTION 

It is well known that Bose-Einstein condensates 
(BECs) with weak attractive interactions between atoms 
can form stable solitons in "cigar-shaped" (nearly one- 
dimensional, ID) traps. In these traps, the gas is strongly 
bound in the transverse plane, while being loosely con- 
fined along the longitudinal axis {z). Using this config- 
uration, stable bright solitons were created in the gas 
of ^Li atoms P, In the condensate of ^^Rb atoms 
trapped in a similar configuration, stronger attraction 
between atoms leads to collapse and emergence of nearly 
3D solitons in a post-collapse state [3]. 

The strongly elongated (cigar-shaped) settings are de- 
scribed by effectively ID equations which were derived, 
by dint of various approximations, from the full 3D 
Gross-Pitaevskii equation (GPE) The deriva- 

tion assumes an ansatz factorizing the 3D wave func- 
tion into a product of a transverse mode and an ax- 
ial (one-dimensional) wave function. As shown in Refs. 
@i @j the substitution of the factorized ansatz in the 
underlying cubic GPE leads, in the general situation, to 
a nonpolynomial Schrodinger equation (NPSE) for the 
axial wave function (or a system of coupled NPSEs for 
a two-component BEC Q). In the case of weak non- 
linearity, the NPSE can be expanded, which leads to a 
simplified ID equation with a combination of cubic and 
quintic terms, the latter one being always attractive; if 
the cubic term is attractive too, the cubic-quintic equa- 
tion gives rise to a family of stable solitons available in 
an exact analytical form Unlike the ordinary ID 

nonlinear Schrodinger (NLS) equation with the attrac- 
tive cubic term, the NPSE, as well as its cubic-quintic 
truncation, may give rise to collapse, which reflects the 
occurrence of the collapse in the underlying cubic GPE 



in three dimensions d, [Tlj . Despite the possibility of 
the collapse, solitons are stable in these models. 

A relevant problem is to derive the NPSE for the cigar- 
shaped trap equipped with a periodic potential, which is 
created in the experiment as an optical lattice (OL), i.e., 
an interferences pattern, by a pair of counterpropagating 
laser beams illuminating the trap in the axial direction 
(the BEC dynamics in periodic potentials was recently 
reviewed in Ref. [ll|). A ID equation of the NPSE 
type including the OL potential was recently derived in 
Ref. [l^. Using that equation, the influence of the peri- 
odic potential on the collapse threshold for axially local- 
ized states in the quasi-lD trap was investigated, both 
for single-peak solitons found in the semi-infinite gap of 
the periodic potential, and multi-peaked solitons found 
in finite bandgaps. The results were compared to direct 
numerical solutions of the full 3D GPE, showing good 
agreement. 

In the experiment, the transverse potential which con- 
fines the atomic gas to the cigar-shaped configuration 
may be axially nonuniform, which corresponds to the cor- 
responding trapping frequency being a function of the 
axial coordinate, ^l±_ — ^^{z) (generally speakin g, i t 
may also depend on time). As proposed in Ref. 
a specially designed nonuniformity {axial modulation) of 
r2_L may be used as an alternative tool for the control 
of dynamics of nearly ID solitons, inducing an effective 
potential for them. In that work, the effective potential 
for a soliton with norm (scaled number of atoms) N was 
found in the limit of weak nonlinearity and long-scale 
axial modulation, 



(l + 5)7Vf]x(C)+0(7'A^'), (1) 



where 7 is an effective nonlinearity coefficient, see Eq. ^ 
below (we will use normalization tantamount to setting 
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TV = 1), C is the coordinate of the soHton's center, and 
integer S is possible intrinsic vorticity of the quasi- ID 
sohton, which is defined below in Eq. ([8]) (ordinary soli- 
tons correspond to 5' = 0). Expression ^ was derived 
by means of a variational approximation (VA) applied 
directly to the full 3D energy functional, cf. Eq. ^ 
below. 

Of special interest is the case of periodic modulation 
of the transverse trapping frequency, with wavenumber 
k, modulation depth a < 1, and amplitude lu±: 



n\{z) ^ Lo]_ [I - acos(2fcz)] . 



(2) 



As suggested by Eq. ([T]), the periodic modulation may re- 
place the OL potential. In particular, this setting may be 
especially relevant to a situation when the quasi- ID mag- 
netic trap is not rectilinear, but rather circular (toroidal). 



which was realized in the experiment [14|. Indeed, the 
OL cannot be created in such a setting, but the induc- 
tion of an effective periodic potential by means of the 
modulation of the transverse trapping frequency is quite 
feasible. 

The objective of the present work is to derive an effec- 
tive NPSE for the quasi-lD trap subject to the periodic 
modulation as per Eq. ^ , and to investigate various self- 
trapped states in this geometry, both delocalized ones 
and solitons. The paper is organized as follows. The 
NPSE is derived, in a general form, in Sec. II. Then, 
the model with the repulsive nonlinearity is considered 
in Sec. III. At first, the analysis includes an additional 
parabolic trapping potential acting in the axial direc- 
tion. Then, this potential is removed, and we analyze 
solutions demonstrating a transition from delocalized so- 
lutions to a gap soliton. The solutions predicted by the 
effective ID equation are compared to their counterparts 
found from the full 3D GPE (bandgap structures gen- 
erated by linearized versions of both equations are com- 
pared too). Section IV is dealing with the case of the 
attractive nonlinearity. The corresponding ground-state 
solitons are found by means of the VA and in a numeri- 
cal form, which are also compared with results produced 
by the full 3D equations. Because the attractive nonlin- 
earity may give rise to collapse, the collapse threshold is 
considered in detail too, by means of both the VA and 
numerical methods. In addition to the ground-state soli- 
tons, which reside in the semi-infinite gap induced by 
periodic modulation ([2]), we also explore solutions with 
the chemical potential belonging to the first finite gap, 
and demonstrate the transition from delocalized states 
to gap solitons in that case. The paper is concluded by 
Sec. V. 



II. THE NONPOLYNOMIAL SCHRODINGER 
EQUATION 

Static BEG configurations can be derived from the 
normalized functional which determines the energy-per- 
atom in terms of order parameter (single-atom wave func- 



tion) "0(1") and includes a generic external axial potential 
V{z): 



-\^^ + 1 [1 _ a cos (2/cz)] {x" + t/2) 



-y(z) + 7r7|0(r)|2]^(r). 



(3) 



Here 7 = 2asN/a± is the adimcnsional strength of the 
interaction between atoms, with the inter-atomic scat- 
tering length and N the number of atoms in the conden- 
sate. In Eq. 13]) lengths are measured in units of the 
characteristic transverse trapping (harmonic-oscillator) 
length, a± = ^Jh/{mu!±), where m is the atomic mass, 
the energy is in units of Tilj^, modulation wave number 
k is in units of a^^, and the wave function is subject to 
the ordinary normalization. 



|V'(r)pdr EE 1. 



(4) 



The chemical potential corresponding to Eqs. ^ and 
O is 



(5) 



Due to the above normalizations, ojx and N are not ex- 
plicitly present in expressions ^ and (O. 

An accurate investigation of the present setting can 
be performed by using the approach which was first de- 
veloped for the GPE with the unmodulated transverse 
trap; after averaging the full 3D equation in the trans- 
verse plane, it leads to an effective one-dimensional NPSE 
[Hi- The derivation of the NPSE starts with the factor- 
ization of the 3D wave function into a product of an arbi- 
trary complex axial wave function, /(z), and the ordinary 
transverse Gaussian ansatz with transverse width 17(2:), 



1 



■ exp 



(x2 



/^a{z) i 2a2(z) 
/(z) being subject to normalization condition 

dz\f{z)\' = l. 



(6) 



(7) 



Gonfigurations including the above-mentioned intrinsic 
vorticity, characterized by positive integer S, correspond 
to the following generalization of Eq. ^ : 



V'(r) 



1 



-p exp 



P 



2a^z) 



tse\f{z), 



(8) 

where p = \/ + y"^ and 9 are the polar coordinates in 
transverse plane {x,y). Nearly-ID solitons with intrinsic 
vorticity were studied, by means of methods similar to 
those used in the present work (although without deriv- 
ing a closed- form NPSE for that case) in Ref. fis']; for a 
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limit case of a strongly elongated 3D trap , similar local- 
ized states were also studied in Ref. [la]. In this work, 
we focus on the ordinary solitons, with = 0. 

Inserting expression ^ in Eq. ([3]), one obtains, ne- 
glecting the z-derivative of a, the following effective (ID) 
energy functional: 



[1 - acos i2kz)] a\z)) + i-^^\fiz)\'^ f{z) . (9) 

The minimization of this functional with respect to f{z), 
taking normalization ^ into regard, leads to the station- 
ary NPSE, 



1 



1 



2 \g'^{z) 



+ [1 - acos (2fcz)] a^{z) 



a\z) 



l/(^)P 



(10) 



where the chemical potential /i appears as the Lagrange 
multiplier generated by the normalization condition. 

An equation for the transverse width is obtained by 
the minimization of fmictional ^ with respect to 17(2;): 



o\z) = 



' l + 7l/(^)P 
1 — a cos {2kz) 



(11) 



The substitution of this expression in Eq. (|T0|) leads to 
a closed-form stationary NPSE for /(z), although in a 
rather cumbersome form. 

In the limit of zero scattering length, 7 = 0, Eqs. (fTO|) 
and (jlip reduce to the linear Schrodinger equation. 



I 52 

+ ^i^) + Vl - acos (2fcz) 



with the effective axial potential. 



Vcff(z) = V{z) + yjl - acos (2fcz) 



(12) 



(13) 



(note that, even in this limit, the full 3D GPE is not 
separable, due to the modulation imposed on the trans- 
verse trapping potential). If the nonlinear term 7|/(z)p 
is small, the NPSE may be approximated by the cubic 
NLS equation with the same effective axial potential and, 
in addition to that, with a periodically- modulated non- 
linearity coefficient: 



1^ 
29^ 



+ K=ff(z) 



+7 v/1 - « COS (2fcz) |/(z)|2j f{z) = ^,f{z). (14) 

Essentially the same result as given by Eq. (HH) was 
obtained, in the simplest approximation, in Ref. [10|, see 
Eqs. (m and (HI). 



III. THE MODEL WITH REPULSIVE 
NONLINEARITY 

A. One-dimensional solutions 

In the case of the repulsive inter-atomic interaction, 
i.e., positive as and 7, the application of the Thomas- 
Fermi (TF) approximation, which neglects the second 
derivative, to Eqs. pI7|) and (fTT|) yields an analytical 
expression for the normalized atomic density, p{z) = 
l/(^)P: 



P(^) 



2 
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A^off(^) - 3 + Aioff(z)y^Aioff(2) + ' 



where an effective local chemical potential is 



tJi-V{z) 



yjl ~ a cos {2kz) 



(15) 



(16) 



Equations P3|) and pB]) generalize the result obtained 
in Ref. Q for the repulsive BEG under the unmod- 
ulated transverse confinement (a = 0). In the limit 
case of strong nonlinearity, 7p 3> 1, Eq. (jlSp reduces 
to p{z) = (4/97) //gg(z), and in the opposite limit of 
7P <C 1, which implies /Xcff (z) — 1 <C 1, it takes the form 

p(z) = (l/7)(Meff(^)-l). 

To obtain accurate results, we solved the time- 
dependent variety of Eq. (|10p (with /i replaced by id/dt 
and, accordingly, d/dz replaced by d/dz), combined with 
Eq. (jlip (without any change in the latter equation) 
numerically, by means of the finite-difference Crank- 
Nicholson method in imaginary time, following the ap- 
proach elaborated in Ref. |13| ■ The explicit axial poten- 
tial was chosen in the usual form, V{z) = (A z) /2, where 
A = LOz/ioj^, and is the axial-confinement frequency. 
In this way, the profiles for p(z), displayed by dashed 
lines in Fig. 1, have been obtained for different values 
of modulation depth a and fixed nonlinearity strength, 
7 = 20. The so obtained NPSE profiles are compared to 
those produced by the TF approximation (dotted lines 
in Fig. 1), see Eq. psp . In the absence of transverse 
modulation, a = (the upper panel in Fig. 1), the axial 
density profile is well approximated by the TF formula. 
At a 7^ (the central and lower panels in Fig. 1), the TF 
approximation much overestimates the density contrast 
between points of local minima and maxima of the total 
axial potential. 



B. Three-dimensional solutions 

It is necessary to compare results produced by the 
NPSE to their counterparts found from a direct numer- 
ical solution of the full GPE in three dimensions. The 
latter equation is obtained by the minimization of energy 
functional ([3]). Figure 1 shows that the density profiles 
generated by the NPSE (dashed lines) are very close to 
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FIG. 1: The axial profile of density p{z) in the self-repulsive 
BEC with 7 = 20, A = 0.1, fc = 1, at different values of the 
modulation amplitude a. Solid lines: results obtained from 
the full three-dimensional GPE, derived by the minimization 
of the underlying energy functional, which is given by Eq. 
(131). Dashed lines: the profiles produced by the numerical 
solution of the one-dimensional nonpolynomial Schrddinger 
equation (NPSE), Eq. (llOp . Dotted lines: results obtained 
from the Thomas-Fermi approximation applied to the NPSE, 
i.e. Eq. (|15p . When using the 3D equation, in this figure and 
below, the axial density displayed in the plots is defined by 
integration in the transverse plane, p(z) = J ^ \ip{v)\^dxdy, 
while in other cases it is simply |/(2;)|'^. 



ones obtained from the 3D equation (solid lines), unless 
a is very large. It is noteworthy that the NPSE gives 
very accurate results for a model with non-separable po- 
tential. 

One may expect that the effective axial periodic po- 
tential induced by the transverse modulation (without 
the inclusion of the axial parabolic trap) may support 
quasi-periodic Bloch states and gap solitons in the self- 
repulsive BEC, as in the case of the ordinary (direct) 
periodic potential [H, [11] . To consider this possibihty, 
a self-consistent numerical method was used, with peri- 
odic boundary conditions. We employed a spatial grid of 
1025 points, covering the interval of —50.26 < z < 50.26, 
which corresponds to 32 periods of the external modu- 
lation. To test the numerical scheme, we have checked 
that the lowest-energy state in the semi-infinite bandgap 
(i.e., the ground state of the system), produced by this 
method, is identical to that found above by the integra- 
tion of the time-dependent NPSE in imaginary time. 

In the upper panel of Fig. 2 we plot, as a func- 
tion of nonlinearity strength 7, the first 50 eigenvalues 
jij, as found by means of the above-mentioned numer- 
ical method from Eqs. HU]) and (HIl) with V{z) = 
and periodic boundary conditions. The lowest 32 eigen- 
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FIG. 2: Upper panel: first 50 eigenvalues m found from the 
NPSE, Eq. Iinil, with V{z) = 0, a = 0.5 and fc = 1, as 
functions of nonlinearity parameter 7. Lower panel: the axial 
density profile, 1/32(2)!^, of the 32th eigenfunction, at a = 0.5, 
for three values of 7. The figure shows that the 32th eigenstate 
becomes a gap soliton at sufficiently large 7. 



states fj{z) belong to the first nonlinear band induced 
by the periodic modulation in the full NPSE, and the 
other states, which are well separated by a bandgap, 
form a second nonlinear band. For some of the intra- 
band states, the numerical method does not converge to 
a single configuration, but rather oscillates between two 
or three configurations with very close eigenvalues. It is 
also possible that the nonlinear model may give rise to 
intraband states which have no counterparts in the linear 
limit. This challenging issue needs special treatment, and 
will be considered elsewhere. As concerns the nonlinear 
eigenstates for which our presently employed numerical 
method leads to convergence (including the gap soliton, 
see below), we have verified, by direct simulations of the 
time-dependent NPSE in real time, that they all are sta- 
ble. 

As a typical example, in the lower panel of Fig. 2 we 
plot the density profile of the 32th state (it has number 
32 in the set of 50 numerically found states) for a — 0.5 
and three different values of 7. One can check that, in 
the linear approximation, this state lies at the top of the 
first Bloch band. With the increase of the nonlinearity 
strength 7, the energy of the 32th state grows, and, in 
doing so, it enters the first finite bandgap (as defined in 
the linear approximation) from below. Figure 2 demon- 
strates that, for 7 = 0.2, this state is still fully delocal- 
ized, being thus similar to a Bloch wave, while for 7 = 0.8 
it becomes localized, with a width much smaller than the 
length of the periodic box. Clearly, this solution may be 
identified as a gap soliton. At 7 = 1.4, the gap soliton 
compresses itself (see Fig. 2) into a still narrower state. 

It is also necessary to check that the bandgap structure 
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in the NPSE with the periodically modulated trapping 
potential in the linear limit (7 = 0) is not, by itself, 
an artifact of the approximation (reduction of the 3D 
equation to the ID limit), but a true feature of the full 
3D model. To this end, we solved the eigenvalue problem 
for the full linear GPE in three dimensions. 



w2 1 5' 



1 2 
+ 2^ 



-a (? cos(2fcz) 



(17) 



where acts on x and y, and, as above, = + y"^. 
The transverse part of solutions to Eq. (fTT]) may be taken 
as an eigenstate of the corresponding 2D harmonic oscil- 
lator, with its quantum numbers and my. Defining 
Gi ~ 2k 'AS the smallest reciprocal lattice vector and 
writing Gn = nGi (n is an integer), we find that the 
solution can be written as 

iP{p,z)= Yl C•™,.™„.G„i^™.,™„(p)e^'«+^"'^ 

nix ^niy.Gn 

(18) 

with coefficients Cm^^my.Gn satisfying a linear eigenvalue 
problem. 



1 2 

{m.j: + my + 1) + - (g + Gn) 



mj:,my,q+G„ 



E.m'^,m'y,G' „ „ 
^m^,my,G„ ^m>_^,rn>^-.q+G' = P"j L-m^ ,m„ ,rj+G„ ^ 

m'^.m'y.G' 

(19) 

where g is a wavenumber in the first Brillouin zone, and 
the matrix A is 

.m'^,m' G' 1 , \ fr 

^m.,m„,G„ ^ ^0^[SG',G„+^+6G',G„-^) X {d,n^,m'^ 



ni'+2 













^my^ 



+ ^(7(S + l)(S + 2)<5™„ 

+\/k + i)( 

my + zj o„j^ 

I (\/« + l)«+2)<^m.,r<+2 



,™;,-2 



+ 6,ny,m'^ (20) 









1 




a- 




+ 2 



+ \/ {rrix + l){mx + 2) ^m^,m^-2) | 



We solved Eq. ((19)) by truncating the sum in expres- 
sion (fT8|) to — 5 < Gn < 5. At first, we also truncated 
the summation to m^ = rUy = (i.e., only the contri- 
bution from the ground state of the transverse potential 
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FIG. 3: Comparison between the dispersion laws as obtained 
from the full 3D equation, Eq. (|19|) . and from the effective ID 
equation, Eq. (|10|) . for k — 1, a = 0.5 and 7 = 0. Left panel: 
Eq. p9p was solved by confining the summation to rrix = 
niy = 0. Right panel: the summation range was extended to 
< rux, rUy < 10. 



was taken into account), which yielded the result shown 
in the left panel of Fig. [31 Then we found more accu- 
rate solutions, by extending the truncated summation to 
< mx,my < 10 (i.e., including the contribution from 
excited transverse states). In that case, Eq. pro- 
duces the dispersion relation shown in the right panel of 
Fig. [3 

These results were compared to the dispersion law as 
found directly from the NPSE, Eq. HU]). We observe that 
the two dispersion laws are identical if the 3D analysis 
is confined to m^ — my — 0, while they are different 
if the effect of states with m^^my > are considered. 
Nevertheless, Fig. [3] shows that the first band, the first 
gap, and half of the second band do not alter essentially, 
if the NPSE is replaced by the full 3D equation, which 
justifies the use of the NPSE approximation for values of 
the chemical potential up to the first half of the second 
band. 



IV. THE MODEL WITH THE ATTRACTIVE 
NONLINEARITY 

A. Numerical results 

The model with the attractive inter-atomic interac- 
tions, i.e. 7 < 0, may be expected to generate bright 
solitons. We analyzed this possibility through the nu- 
merical solution of the NPSE equation, in the absence of 
the direct axial potential, V{z) = 0. In Fig. 4, we plot 
axial density profiles of the thus found solitons for differ- 
ent values of modulation depth a and fixed nonlinearity 
strength, g = —7. 
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FIG. 4: Axial density profiles, p{z), of the ground-state bright 
soliton in the attractive model with g = 0.5 and fc = 1, for 
different values of modulation amplitude a. Dashed lines: 
results obtained from the NPSE, Eq. (fTIJ)l . Solid lines: results 
obtained from the 3D equation, derived by the minimization 
of energy functional ^ . 



For a = 0, the soliton's profile, with a single maxi- 
mum, may be well fitted by f{z) = {y/g/2) sech{gz/2), 
which is an asymptotically exact solution (the usual NLS 
soliton) in the above-mentioned weakly nonlinear limit 
corresponding to gf^ ^ 1, provided that z varies on the 
entire real axis. On the other hand, if z belongs to a 
finite interval, — L/2 < z < L/2, with periodic bound- 
ary conditions, f{z + L) — f{z), it is known [l9j| that 
the NPSE with a = yields a spatially uniform ground 
state, f{z) = for sufficiently weak nonlinearity, 

< g < TT^ / L; the ground state develops a spatial struc- 
ture at g > TT^/L . As shown in Fig. 4, for nonzero a 
the soliton profile features several local maxima and min- 
ima due to the action of the effective periodic potential. 
Thus, under such conditions, the Bose condensate self- 
traps into a multi-peaked soliton, which occupies several 
cells of the periodic modulation. 

We also compared the results yielded by the NPSE 
with those found from the direct numerical solution of 
the full GPE in three dimensions. Figure 4 shows that 
the density profiles generated by the NPSE (dashed lines) 
practically coincide with the ones obtained from the 3D 
GPE (solid lines), unless a is very large. 

It is also interesting to analyze the behavior of the 
density profile of the bright soliton as a function of the 
wavenumber k. In Fig. 5 we plot the density profile p{z) 
at g = 0.5 and a — 0.5, for four values of k. The figure 
shows that, at k — 0.5, the profile is strongly localized 
within one single site of the effective periodic axial po- 
tential: accordingly, we call the corresponding solution 
a single-site soliton (cf. Ref. [l^l, where the NPSE for 
the model with attraction and an explicit periodic axial 
potential was considered). At k < 0.75, delocalization 



FIG. 5: Axial density profiles, p{z), of the ground-state bright 
soliton in the attractive model with g = 0.5 and a = 0.5, for 
different values of wave number k. Results obtained from the 
NPSE, Eq. dlD}. 



of the bright soliton, which occupies more than one site, 
is observed. We call it a multi-site soliton (the distinc- 
tion between the strongly and weakly localized solitons 
is also observed in the ID cubic NLS equation with the 
self- attractive nonlinearity (23|). 



B. Variational vs numerical results 

From the numerical results presented above we see that 
the profiles of the localized solutions of Eq. both in 
the single-site (strong attraction, large values of g) and 
in the multi-site soliton cases but with a weak transverse 
modulation (small values of a) are smooth. This obser- 
vation suggests that one could achieve some analytical 
insight into the model with attraction by using a VA with 
the Gaussian ansatz (a review of the VA can be found in 
Ref. 0), 



^(r) 



1 



7r3/4o-7;i/ 



YT^exp 



2cr2 



exp 



(21) 

where a and 77 are variational parameters accounting for 
the transverse and axial width of the BEG. Inserting this 
ansatz in Eqs. ^ and ((5]), with V{z) = (and 7 = —g), 
we obtain the respective expressions for the energy-per- 
atom functional and chemical potential, 

^ = i+4^ + ^ [1 - aexp(-fcV)] -^-^(^^ - 

(22) 

(23) 



/I = 



9 
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Next, we minimize the energy, demanding dE/da 
dE/drj = 0, which yields the variational equations, 
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[l - aexp(-fc2r?2)] a"^ = I - 



(27r)i/2,y 



(271)1/2^4 



(24) 



(25) 



which can be easily solved numerically [2l| . The solutions 
provide for a minimum of the energy only if the curvature 
of the energy surface, i?(ry, ti), is positive, i.e.. 



d^Ed^E / d^E 
drf 9(7^ \drida 



> 



(26) 



As concerns the dynamical stability of the solitons 
against small perturbations, it may be, first of all, es- 
timated by means of the VK criterion [22|. According 
to it, a necessary stability condition is dfi/dg < (in the 
present notation) , if the soliton family is described by de- 
pendence /i(g) (note that the VK criterion does not apply 
to gap solitons in the model with the repulsive nonlinear- 
ity, therefore it was not used in the previous section). 

In Fig. 6, we plot axial length rj and transverse width a 
versus interaction strength g for k = 1 and four different 
fixed values of modulation parameter a. The figure shows 
that the soliton in the self-attractive BEC exists up to a 
critical value of the nonlinearity strength, gc- At g > gc, 
the 3D collapse of the nearly-lD soliton occurs, which 
is a well-known result in the case of a — 0, S 0- 
Note that, for a = 0, the VA predicts gc — 1.55, which 
is somewhat higher than gc = 1.34 obtained from the 
numerical solution of the full GPE in three dimensions 

The axial length of the soliton, rj, diverges as g drops 
to zero, while the transverse width a approaches 1, ac- 
tually becoming equal to the above-mentioned harmonic- 
oscillator length, a±. On the other hand, as g approaches 
gc, both r] and a remain finite and smaller than 1. 

New results are presented in Fig. 6 for a (recall 
previous works were only dealing with the case of a = 0). 
At small a, the figure shows only a slight distortion of 
the curves. A qualitative change is observed at a > 0.36, 
when there appear two stable branches for both a{g) and 
r]{g), the curves for r]{g) displaying a clear gap. The lower 
branches of the r]{g) and cr(g) dependences exists only in 
a finite interval, which we denote as 



gni< g < gc, 



(27) 



while the upper branches extend up to 5 = 0, in inter- 
val < g < gu, with gM < gc- Physically, the lower 
branches (with smaller values of axial length rj) corre- 
spond to single-site solutions, where the self-attractive 
BEC is strongly localized - essentially, within a single 
cell of the modulation structure. The upper branches 
of 77(5) and cr{g) correspond instead to the weakly local- 
ized solutions (with a larger axial length), which occupy 
several cells (sites). 
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FIG. 6: Transverse width a (dashed lines) and axial length 
rj (solid lines) of the quasi- ID bright soliton in the model 
with attraction, versus the nonlinearity strength, g = —7 = 
2\as\N/a±, as predicted by the variational approximation 
based on ansatz ((21]) and Eqs. ((24)) . ((25|) . On the right side, 
all curves terminate at the critical point, g = gc, past which 
the collapse occurs. The dependences are displayed for k = 1 
and different values of the modulation depth a. 



Analysis of expressions ((H|) and ((53|) demonstrates 
that, in interval ((77[) . the multi-site and single-site so- 
lution may assume the role of ground state (the one cor- 
responding to the lowest energy). However, this analysis 
also suggests that both families are dynamically stable, 
as they always meet the VK criterion, d^/dg < 0. Direct 
numerical simulations (not shown in detail here) have 
confirmed this conjecture. 

For a > 0.86, the numerical solution of Eqs. ((M(l and 
((^5]) demonstrates that the lower border of existence in- 
terval for the single-site soliton, g,„, vanishes, but 
this is an artifact of the VA, which occurs in other con- 
texts too ^2§\. It is explained by the above-mentioned 
inadequacy of the Gaussian ansatz in the limit of weak 
nonlinearity, i.e., for widely spread small-amplitude soli- 
tons featuring a multi-peaked shape. In this situation, 
one may, in principle, apply a more sophisticated ansatz, 
combining the Gaussian and periodic functions, such as 
cos(2fciz); however, the generalized ansatz results in a 
cumbersome algebra [2^ , therefore we do not pursue such 
an approach here. 

In the upper panel of Fig. 7, we plot the density pro- 
file, p(z), of the soliton for different values of the self- 
attraction strength, g, and fixed modulation parameters, 
a = 0.9 and k = 1. As seen in the figure, the increase 
of g may strongly compress the soliton in the axial direc- 
tion, making the secondary maxima very small. In this 
case, the condensate actually self-traps into a single-peak 
soliton, which occupies only one cell of the modulation 
structure. 

Contrary to the numerical solution of the NPSE which 
shows a smooth crossover between multi-site and single- 
site solitons, the VA predicts a discontinuous transition. 
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FIG. 7: Upper panel: the axial density profile, p(z) — \f{z)\ , 
of the soliton in the attractive model, with a = 0.9 and fe = 1, 
obtained for different values of the self-attraction strength g 
from the numerical solution of the NPSE, Eq. (|10|l . Lower 

panel: axial length of the soliton, (^^)^ ^, as a function of p, 
for a — 0.9 and k — 1. The lower panel includes results pro- 
vided by the variational approximation based on the Gaussian 
ansatz, see Eqs. (|24|l and (|25|l . and those obtained from the 
numerical solution of the NPSE. 



In the lower panel of Fig. 7, we plot the axial length, 
•y/ (z^), of the ground-state bright soliton as a function of 
strength g, for a = 0.9 and k = 1. The jump predicted by 
the variational calculation at g ~ 0.5, is a consequence of 
the inadequacy of of ansatz (|2ip . which assumes the sim- 
ple Gaussian waveform for the axial wave function, to de- 
scribe multi-peaked states, as was also recently shown in 
the study of the quasi- ID model with the self- attraction 
and axial optical lattice p^ . 

As it is well known, cold Bose atoms with attractive 
interactions collapse if the interaction strength exceeds a 
threshold value, gc- Obviously, as g increases towards gc, 
the profile becomes narrower and narrower, and therefore 
close to collapse in the presence of a transverse modula- 
tion, its shape should correspond to a single-site soliton, 
where the VA is appropriate. It is thus interesting to 
predict the collapse threshold, gc, as a function of pa- 
rameters a and k of the transverse modulation by using 
our gaussian variational ansatz. In Fig. 8, we display the 
dependence gdk) predicted by the VA at five fixed val- 
ues of modulation depth a. For given a, the critical value 
gc has its maximum at fc = 0, which is natural, as Eq. 
([2]) yields, in this case, the smallest constant value of the 
transverse-trapping frequency. Equation Q also helps 
to understand another feature observed in Fig. 8., viz., 
that gc slowly diverges at fc = as a approaches 1 (for 
instance, gc = 4.88 at a = 0.99). Note also that there 
exists a modulation wavenumber, kc, at which gc attains 
its minimum, i.e., the collapse has the lowest threshold. 
Figure 8 shows that this minimum decreases with the in- 
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FIG. 8: Critical strength Qc for the collapse of the quasi-lD 
soliton in the model with attraction versus modulation wave 
number k, as predicted by the variational approximation. The 
dependences are displayed at fixed values of modulation depth 
a. 



crease of a, which may be understood too: as mentioned 
above, the strong potential tends to squeeze the entire 
condensate into a single cell of the modulation structure, 
which facilitates the onset of the collapse. On the other 
hand, at large values of fc, gc becomes asymptotically con- 
stant, as the interaction of the condensate with the short- 
period modulation becomes exponentially weak, hence it 
produces little effect on the collapse threshold. 

In Table 1, we display the comparison of the VA- 
predicted critical value, gc, to results following from the 
numerical solution of NPSE pO)) (for a = 0.9). The vari- 
ational and numerically found critical values are seen to 
be in qualitative agreement. Note that, at fc = 0, the crit- 



ical value for the NPSE is gc 
is an obvious extension of gc - 



-1/4 



which 



.(4/3)(l-a) 
4/3 found before at a = 



fc 


9c 


(var ) 

9 c 




a(0) 


0.00 


2.12 


2.76 


2.16 


0.86 


0.10 


1.99 


2.52 


1.20 


0.72 


0.25 


1.74 


2.17 


0.91 


0.65 


0.50 


1.55 


1.85 


0.75 


0.56 


1.00 


1.41 


1.57 


0.76 


0.57 


1.50 


1.34 


1.47 


1.19 


0.79 


2.00 


1.32 


1.48 


1.25 


0.83 


2.50 


1.45 


1.54 


0.99 


0.75 


3.00 


1.44 


1.55 


0.93 


0.73 


3.50 


1.42 


1.55 


0.92 


0.73 



Table 1. Properties of the bright soliton in the model with 
attraction in a proximity to the collapse, as found from the 
numerical solution of Eq. (|1U|I at a = 0.9 and different val- 
ues of modulation wavenumber k: is the critical value of 
the nonlinearity coefficient at the collapse point, \J < z'^ > is 



the cixial width of the soUton, and (j(0) its transverse width 
(at 2 = 0) at the same critical point. For comparison, also 
included are values of Qa predicted by the variational approx- 

,. (var) 

imation, gc 

Table 1 demonstrates that, with the increase of fc, the 
critical value, Qc, drops from the largest value, corre- 
sponding to fc = 0, to a minimum at fc = 2, and then 
slightly increases with the further increase of k (the mini- 
mum at fc = 2 seems quite shallow from the side of fc > 2). 
This feature and the related ones can be explained if one 
notices that fc = corresponds to a constant value of the 
modulation factor in Eqs. PH)) and pT|) . 



/3fc=o = - acos (2fcz)|fc=o = Vl - a, 



(28) 



and, on the other hand, for large fc (formally, for fc — > oo), 
the modulation factor should be replaced by its average. 



/3fc= 



E (Vl-acos(2fcz)) = ^VTT^E [^|^^ , 

(29) 

where E is the complete elliptic integral of the second 
kind. Then, Eqs. PH)) and (fTT|) with the constant modu- 
lation factor (3 (and without the extra potential, = 0), 
take the form of 




1 



x/l-5l/WP 



^ (3 ^ 



(30) 



(recall g = —7). Further, the coefficient /3 can be elim- 
inated from Eq. (|30|) by means of an obvious rescaling 
[which also takes into regard the condition that the nor- 
malization of the solution must keep the form of Eq. ^] : 

As a consequence, the critical values of g, together with 
the respective length scale, which correspond to the dif- 
ferent constant values of (3, are related as follows: 



f3k=oc 



(31) 



For a = 0.9 (the value for which numerical data are col- 
lected in Table 1), Eqs. ^ and ^ yield Po ~ 0.32 
and Poo ~ 0.93, hence the ratios predicted by Eq. (^1]) 
take the value ~ 0.59. On the other hand, the numeri- 
cal data from Table 1 yield {gc)k=2 I i9c)k=o ~ 0.62, and 
( -v/Tz^y ] / ( y/ (z^) ] « 0.58. Thus, approximating 

V / fc=2 V / k=0 

(7c at fc = 2 by {gc)k^ao (recall the change of at fc > 2 
is insignificant), one can explain (at least, in a crude ap- 
proximation) effects caused by the transition from long- 
scale to short-scale spatial modulation. 
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FIG. 9: First 50 eigenvalues found from the numerical 
solution of NPSE, Eq. ([TOj, with a = 0.5 and fc = 1, as 
functions of the nonlinearity strength in the model with at- 
traction, g = —7. 



C. Solitons in the first finite bandgap 

The above analysis was dealing with solitons (in the 
attractive model) whose chemical potential belongs to 
the semi-infinite bandgap in the linear spectrum of Eq. 
(fTU)) . On the other hand, it is known that the cubic 
self-attractive nonlinearity may also give rise to solitons 
located in higher-order (finite) bandgaps. Here, we re- 
port soliton solutions of the latter type, found from the 
NPSE by means of the self-consistent numerical method 
with periodic boundary conditions, which was outlined 
in the previous section. 

In Fig. 9 we plot the first 50 eigenvalues fij, as found 
from the numerical solution of the stationary nonlinear 
NPSE, versus the nonlinearity strength g. The first 32 
eigenstates form the first band, and the other 18, which 
are well separated by a gap, cluster into the second band. 
Figure 9 shows that the lowest eigenvalue and the 33rd 
one split off from the first and second bands and move 
down (till the onset of the collapse) with the increase of 
g, thus giving rise to localized states in the semi-infinite 
and first finite gaps, respectively. It is worthy to note 
that the second eigenvalue, which originally belonged to 
the first band, also splits off from it at larger values of g. 
We have verified that the corresponding nonlinear eigen- 
state become localized, as one may expect. Qualitatively 
similar findings were reported in Ref. [27j , which was 
dealing with a numerical solution of the ordinary cubic 
GPE in one dimension, and, more recently, in Ref. [T3| 
which was dealing with the self-attractive BEG in the 
quasi-lD trap with an axial OL potential. 

The density profiles, 1/33(2)^, of the 33th state which 
develops into a soliton belonging to the first finite 
bandgap, are displayed, for a = 0.5 and six different val- 
ues of interaction strength g, in Fig. 10. The figure shows 
that this state is still fully delocalized (being similar to 
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FIG. 10: Axial density profiles, | /as (2;)!^, of the 33th nonlinear 
eigenfunction of the NPSE (which develops from a Bloch-like 
delocalized state into a gap soliton) for a = 0.5, k = 1 and 
six values of g. 



a Bloch wave) at g = 0.4, while at 5 = 0.6 it becomes 

localized, with the width much smaller than the length 
of the periodic box, featuring many local maxima and 
minima (zeros). With further increase of g, the gap soli- 
ton keeps compressing itself. Comparing Figs. 4 and 10, 
we conclude that the ground-state solitons, which reside 
in the semi-infinite gap, are drastically different from the 
gap solitons, i.e., ones found in the first finite bandgap. 
First, for the same parameters, the ground-state solitons 
are localized much stronger, and their local density min- 
ima are not zeros, unlike those of the gap solitons. In 
addition, it is noteworthy that local density maxima of 
the groimd-state solitons correspond to minima of the ef- 
fective periodic potential, while the local maxima of the 
gap soliton correspond to maxima of the periodic poten- 
tial. 

Although the gap solitons cannot play the role of 
ground states, they, as well as the ground-state solitons, 
are stable in direct simulations of the time-dependent 
variant of the NPSE equation. Therefore, the gap-soliton 
states are relevant to the experiment. 



V. CONCLUSIONS 

In this work, we have derived the effective one- 
dimensional NPSE (nonpolynomial Schrodinger equa- 



tion) for a cigar-shaped trap whose transverse confining 
frequency is periodically modulated along the axial direc- 
tion, thus inducing an effective periodic axial potential. 
Besides the usual quasi- ID linear geometry, the model 
may also be relevant as a means of creating an effec- 
tive periodic potential in toroidal traps. In both cases 
of the repulsive and attractive nonlinearity, delocalized 
states and solitons were found, by means of numerical 
methods (which were applied to both the effective NPSE 
and the underlying 3D Gross-Pitaevskii equation) and 
VA (variational approximation; this method was applied 
to ground-state solitons in the model with attraction, and 
to the prediction of the collapse threshold in this case). 
It was found that the numerical solution to the NPSE is 
always extremely close to the full 3D solutions. The VA 
yields quite reasonable results too, except for the descrip- 
tion of the crossover from single-site to multi-site solitons: 
numerical results reveal that the crossover is smooth and 
does not include a jump, contrary to the prediction of the 
VA. This shortcoming of the VA is explained by the fact 
that the simple Gaussian ansatz, on which the approxi- 
mation is based, cannot adequately grasp the transition 
that alters the shape of the soliton, giving it the multi- 
peaked structure. The transition from delocalized states 
to gap solitons was studied in detail (by means of numer- 
ical methods) in the first finite bandgap, for both cases 
of the repulsive and attractive nonlinearities. 

The above results, presented in terms of the dimen- 
sionless equations, can be easily translated into physical 
units. For instance, by considering an attractive Bose- 
Einstein condensate made of ^Li atoms, with scattering 
length Us = —1.45 nm, and choosing the transverse con- 
fining frequency as uj± =2Tr x 100 Hz, we have a typical 
value of the modulation period, A = na^/k = 12 /im, for 
k = 1 (recall it was a typical value for which the results 
were reported above). In this case the critical number of 
atoms at the collapse threshold is Nc = gcO-±/{2\as\) ~ 
1850. 

The results may be quite useful to design new exper- 
iments in toroidal traps, where an azimuthal periodic 
potential can be induced, as explained above, by the 
spatially-modulated transverse confinement. 
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